double a1 = 10e-6;
double b1 = 1500e-6;

double a2 = 10e-6;
double b2 = 1500e-6;

double ft(double *x, double *par) {

  double t = x[0];
  double tau_H = par[0];
  double tau_G = par[1];
  double N_H = par[2];
  double N_G = par[3];

  double binW = par[5];

  double pH = N_H * TMath::Exp(-t / tau_H) / ((tau_H * (TMath::Exp(-a2 / tau_H) - TMath::Exp(-b2 / tau_H))));
  double pG = N_G * TMath::Exp(-t / tau_G) / ((tau_G * (TMath::Exp(-a1 / tau_G) - TMath::Exp(-b1 / tau_G))));

  return binW * (pH + pG) + par[4];
}

double ft_Gd(double *x, double *par) {

  double t = x[0];
  double tau_G = par[0];
  double N_G = par[1];
  double binW = par[2];

  double pG = N_G * TMath::Exp(-t / tau_G) / ((tau_G * (TMath::Exp(-a1 / tau_G) - TMath::Exp(-b1 / tau_G))));

  return binW * pG + par[4];
}

int *h2dGetXY(int k) {
  int *res = new int[2];
  res[1] = 5 - (k - 1) / 5;
  res[0] = ((k % 5) == 0) ? 5 : (k % 5);
  return res;
}

void AddConstant(double val, double err, double &val_base, double &err_base);

void SaveConstant1D_lijj(double X[3][4][5][5], double Y[3][4][5][5], double Y_err[3][4][5][5]);

int Fit_time_basic() {

  gStyle->SetPaintTextFormat(".1f");
  TF1 *func_pol = new TF1("func_pol", ft, 0, 1.5e-3, 6);
  func_pol->SetNpx(1e6);
  func_pol->SetLineColor(2);
  func_pol->SetParNames("#tau_{LS}", "#tau_{Gd}", "N_{LS}", "N_{Gd}", "p0", "binW");
  func_pol->SetLineWidth(2);

  func_pol->SetParLimits(2, 0, 1e8);
  func_pol->SetParLimits(3, 0, 1e8);

  func_pol->FixParameter(0, 210e-6);
  func_pol->FixParameter(1, 30e-6);

  double EdgesX[6] = { 0, 1e6, 2e6, 2.5e6, 3.25e6, 4e6 };
  double EdgesY[6] = { -2000, -1400, -500, 500, 1400, 2000 };

  TString top_dir = "/dybfs2/users/lijj/DATA/Accsub/IBD/IBD_version0.8/phi7000mm_Im_withGridTest/Output/";
  // TString top_dir = "phi800mm_grid_time_neutron/Output/";

  bool flag_Gd = false;

  int nAD[3] = { 2, 2, 4 };

  // fitted constant
  TH2D *h[3][4];

  TH2D *h_N_Gd[3][4];
  TH2D *h_N_LS[3][4];

  double val_constant_bulk[8] = { 0 };
  double err_constant_bulk[8] = { 0 };
  double val_NGd_bulk[8] = { 0 };
  double err_NGd_bulk[8] = { 0 };
  double val_NLS_bulk[8] = { 0 };
  double err_NLS_bulk[8] = { 0 };
  double iad_num[8];
  int i_clock_bulk = 0;

  double constant_X[3][4][5][5];
  double constant_Y[3][4][5][5];
  double constant_Y_err[3][4][5][5];
  double offset_t[3][4] = { { -0.12, -0.09, 0, 0 }, { -0.06, -0.03, 0, 0 }, { 0.03, 0.06, 0.09, 0.12 } };

  for (int i = 0; i < 3; i++) {
    TFile *f = new TFile(top_dir + TString::Format("EH%d/EH%d_tot.root", i + 1, i + 1), "READ");
    for (int j = 0; j < nAD[i]; j++) {
      h[i][j] = new TH2D(Form("h_%d_%d", i + 1, j + 1), "", 5, EdgesX, 5, EdgesY);
      h_N_Gd[i][j] = new TH2D(Form("h_N_Gd_%d_%d", i + 1, j + 1), "", 5, EdgesX, 5, EdgesY);
      h_N_LS[i][j] = new TH2D(Form("h_N_LS_%d_%d", i + 1, j + 1), "", 5, EdgesX, 5, EdgesY);
      h[i][j]->SetXTitle("X^{2}+Y^{2} [mm^{2}]");
      h[i][j]->SetYTitle("Z [mm]");
      h_N_Gd[i][j]->SetXTitle("X^{2}+Y^{2} [mm^{2}]");
      h_N_Gd[i][j]->SetYTitle("Z [mm]");
      h_N_LS[i][j]->SetXTitle("X^{2}+Y^{2} [mm^{2}]");
      h_N_LS[i][j]->SetYTitle("Z [mm]");
      TH1D *h_voxel[25];

      for (int k = 0; k < 25; k++) {
        h_voxel[k] = (TH1D *)f->Get(Form("h1dDt_%d_grid_%02d", j + 1, k + 1));

        double binW = 1e-6;
        if ((k + 1) == 6 || (k + 1) == 7 || (k + 1) == 8 || (k + 1) == 11 || (k + 1) == 12 || (k + 1) == 13 || (k + 1) == 16 || (k + 1) == 17 || (k + 1) == 18) {
          flag_Gd = true;
          func_pol->ReleaseParameter(5);
          func_pol->FixParameter(5, binW);
          func_pol->SetParameters(210e-6, 30e-6, 10, h_voxel[k]->GetBinContent(20), 1, binW);
        } else {
          func_pol->ReleaseParameter(5);
          func_pol->FixParameter(5, binW);
          h_voxel[k]->Rebin(10);
          binW *= 10;
          func_pol->SetParameters(210e-6, 30e-6, h_voxel[k]->GetBinContent(2), 10, 1, binW);
        }

        h_voxel[k]->Fit("func_pol", "NOSQ", "", flag_Gd ? 10e-6 : 10e-6, 1.5e-3);
        TFitResultPtr fptr = h_voxel[k]->Fit("func_pol", "NOSQ", "", flag_Gd ? 10e-6 : 10e-6, 1.5e-3);
        if (((i + 1) == 1 && (j + 1) == 1 && (k + 1) == 14)) {
          //TFitResultPtr fptr = h_voxel[k]->Fit("func_pol", "SQ", "", flag_Gd ? 10e-6 : 10e-6, 1.5e-3);

          //new TCanvas(Form("EH%d-AD%d-voxel_%d",i+1,j+1,k+1),Form("EH%d-AD%d-voxel_%d",i+1,j+1,k+1));
          //h_voxel[k] -> Draw();
          //cout<<Form("voxel %d, C/Ndf= %.2f",k+1,fptr->Chi2()/fptr->Ndf())<<endl;
        }

        int iX = h2dGetXY(k + 1)[0];
        int iY = h2dGetXY(k + 1)[1];

        double N_LS = func_pol->GetParameter(2);
        double N_Gd = func_pol->GetParameter(3);
        double constant = func_pol->GetParameter(4);

        double N_LS_err = func_pol->GetParError(2);
        double N_Gd_err = func_pol->GetParError(3);
        double constant_err = func_pol->GetParError(4);

        constant_X[i][j][5-iY][iX-1]=k+1+offset_t[i][j];
        constant_Y[i][j][5-iY][iX-1]=constant;
        constant_Y_err[i][j][5-iY][iX-1]=constant_err;

        //some invalid result
        if (N_Gd < 1 || N_Gd < N_Gd_err)
          N_Gd = 0;
        if (N_LS < 1 || N_LS < N_LS_err)
          N_LS = 0;

        // Add each voxel into bulk
        AddConstant(constant, constant_err, val_constant_bulk[i_clock_bulk], err_constant_bulk[i_clock_bulk]);
        AddConstant(N_Gd, N_Gd_err, val_NGd_bulk[i_clock_bulk], err_NGd_bulk[i_clock_bulk]);
        AddConstant(N_LS, N_LS_err, val_NLS_bulk[i_clock_bulk], err_NLS_bulk[i_clock_bulk]);

        h[i][j]->SetBinContent(iX, iY, constant);
        h[i][j]->SetBinError(iX, iY, constant_err);

        h_N_Gd[i][j]->SetBinContent(iX, iY, N_Gd);
        h_N_Gd[i][j]->SetBinError(iX, iY, N_Gd_err);

        h_N_LS[i][j]->SetBinContent(iX, iY, N_LS);
        h_N_LS[i][j]->SetBinError(iX, iY, N_LS_err);

        flag_Gd = false;
      }

      // cout<<Form("EH%d-AD%d: %.4f +/- %.4f",i+1,j+1,sum_constant,err_constant)<<endl;
      iad_num[i_clock_bulk] = i_clock_bulk + 1;
      i_clock_bulk++;

      new TCanvas(Form("EH%dAD%d-Constant", i + 1, j + 1), Form("EH%dAD%d-Constant", i + 1, j + 1));
      h[i][j]->Draw("COLZ+TEXT+E");
      h[i][j]->SetMarkerSize(2);
      gPad->SaveAs(Form("Figures/EH%dAD%d-Constant.pdf", i + 1, j + 1));

      new TCanvas(Form("EH%dAD%d-N-Gd", i + 1, j + 1), Form("EH%dAD%d-N-Gd", i + 1, j + 1));
      h_N_Gd[i][j]->Draw("COLZ+TEXT+E");
      h_N_Gd[i][j]->SetMarkerSize(2);
      gPad->SaveAs(Form("Figures/EH%dAD%d-N-Gd.pdf", i + 1, j + 1));

      new TCanvas(Form("EH%dAD%d-N-LS", i + 1, j + 1), Form("EH%dAD%d-N-LS", i + 1, j + 1));
      h_N_LS[i][j]->Draw("COLZ+TEXT+E");
      h_N_LS[i][j]->SetMarkerSize(2);
      gPad->SaveAs(Form("Figures/EH%dAD%d-N-LS.pdf", i + 1, j + 1));
    }
  }

  TFile *g_out_file = new TFile("Fitted_voxel_NNC.root", "RECREATE");
  TGraphErrors *gConstant = new TGraphErrors(8, iad_num, val_constant_bulk, 0, err_constant_bulk);
  TGraphErrors *gNGd = new TGraphErrors(8, iad_num, val_NGd_bulk, 0, err_NGd_bulk);
  TGraphErrors *gNLS = new TGraphErrors(8, iad_num, val_NLS_bulk, 0, err_NLS_bulk);
  gConstant->Write("Constant");
  gNGd->Write("N_Gd");
  gNLS->Write("N_LS");
  g_out_file->Close();

  SaveConstant1D_lijj(constant_X,constant_Y,constant_Y_err);

  return 0;
}

void AddConstant(double val, double err, double &val_base, double &err_base) {
  val_base += val;
  err_base = sqrt(pow(err_base, 2) + pow(err, 2));
}

void SaveConstant1D_lijj(double X[3][4][5][5], double Y[3][4][5][5], double Y_err[3][4][5][5]) {
  //color and offset of 8 ADs
  int color_t[3][4] = { { 1, 2, 0, 0 }, { 3, 4, 0, 0 }, { 9, 6, 7, 8 } };

  TLegend *lg = new TLegend(0.25, 0.45, 0.75, 0.95);
  lg->SetNColumns(2);
  lg->SetTextSize(0.15);
  lg->SetFillColor(0);
  lg->SetFillStyle(0);

  TCanvas *c1 = new TCanvas("c1", "c1", 800, 1000);
  c1->Divide(1, 6, 0.0, 0.0);

  int nAD[3] = { 2, 2, 4 };

  TGraphErrors* gr_Mean[3][4][5];
  
  for (int isite = 1; isite < 4; isite++) {
    for (int AdNo = 1; AdNo < nAD[isite - 1] + 1; AdNo++) {

      for (int i = 0; i < 5; i++) {
        c1->cd(i + 1);
        gr_Mean[isite - 1][AdNo - 1][i] = new TGraphErrors(5, X[isite-1][AdNo-1][i], Y[isite-1][AdNo-1][i], 0, Y_err[isite-1][AdNo-1][i]);
        gr_Mean[isite - 1][AdNo - 1][i]->Draw((isite == 1 && AdNo == 1) ? "AP" : "P");
        gr_Mean[isite - 1][AdNo - 1][i]->SetMarkerColor(color_t[isite - 1][AdNo - 1]);
        gr_Mean[isite - 1][AdNo - 1][i]->SetLineColor(color_t[isite - 1][AdNo - 1]);
        gr_Mean[isite - 1][AdNo - 1][i]->GetXaxis()->SetTitle("voxel number");
        gr_Mean[isite - 1][AdNo - 1][i]->GetYaxis()->SetTitle("Relative difference [%]");
        gr_Mean[isite - 1][AdNo - 1][i]->GetYaxis()->SetNdivisions(508);
        gr_Mean[isite - 1][AdNo - 1][i]->GetYaxis()->SetRangeUser(-50, 50);

        if (i == 0)
          lg->AddEntry(gr_Mean[isite - 1][AdNo - 1][i], Form("EH%d-AD%d", isite, AdNo), "lpe");

        gPad->SetGridx();
        gPad->SetGridy();
        gPad->SetBottomMargin(0.2);
        gPad->SetTopMargin(0.0);
        if (i == 4)
          gPad->SetBottomMargin(0.2);
        gr_Mean[isite - 1][AdNo - 1][i]->GetXaxis()->SetLabelSize(0.14);
        gr_Mean[isite - 1][AdNo - 1][i]->GetYaxis()->SetLabelSize(0.14);
        gr_Mean[isite - 1][AdNo - 1][i]->GetXaxis()->SetTitleSize(0.001);
        gr_Mean[isite - 1][AdNo - 1][i]->GetYaxis()->SetTitleSize(0.001);
      }
      //
    }
  }
  c1->cd(6);
  lg->Draw("same");
  c1->cd();
  c1->SaveAs("Figures/summary_constant.eps");
}
